Supersymmetric Langevin equation to explore free energy landscapes 
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The recently discovered supersymmetric generalizations of Langevin dynamics and Kramers equa- 
tion can be utilized for the exploration of free energy landscapes of systems whose large time-scale 
separation hampers the usefulness of standard molecular dynamics techniques. The first realistic 
application is here presented. The system chosen is a minimalist model for a short alanine peptide 
exhibiting a helix-coil transition. 
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The problem of identifying and exploring critical re- 
gions on free energy landscapes is central to many disci- 
plines lying at the interface between physics and chem- 
istry [l| : It is an outstanding issue in the study of atomic 
clusters, glasses, supercooled liquids, biopolymer dynam- 
ics, and protein folding. 

Many authors 0, S, i, i, 0, 0, H have been develop- 
ing methods to find reaction paths based on a statistical 
description of the ensemble of pathways connecting cer- 
tain phase space regions. A radic ally different approach 
has been recently proposed [l^, [ij : The supersymme- 
try hidden in the Kramers equation can be made explicit 
by coupling to the "bosonic" degrees of freedom in the 
phase space an equal number of "fermionic" ones in the 
tangent space. As the usual Kramers equation is the base 
of molecular dynamics (henceforth simply MD) simula- 
tions, so its supersymmetric generalization gives rise to 
an enhanced MD (let's call it SuSy MD), suitable for the 
study of systems characterized by time-scale separation. 
Rigorous theoretical arguments and 1- and 2-dimensional 
toy model applications have established that, in a purely 
energetic landscape, SuSy MD is able to find reaction 
paths in a time much shorter than the activation time. 

However, in order to pave the way to realistic problem 
applications, the feasibility of SuSy MD must be estab- 
lished in the framework of high dimensional free energy 
landscapes. By studying a model not very computation- 
ally demanding, but possessed of all the other complica- 
tions which characterize current research in biomolecular 
simulations, we provide the missing piece of evidence that 
the new method is able to signal the presence of entropic 
as well as energetic barriers. We show that the iden- 
tification of the reaction path obtained by SuSy MD is 
fully consistent with the results of standard MD, with 
the advantage that the simulation time needed is orders 
of magnitude shorter. 

In the following we first briefly review the super- 
symmetric formulation of Kramers equation, urging the 
reader interested in more details to peruse Refs. [lCll. [T]|: 
With respect to the existing literature, however, we try 
to clarify the peculiar problems raised by the application 
to free energy landscapes. A brief description of techni- 
cal aspects of the method used is followed by the results 
of its application to our test system. 



I. SUPERSYMMETRIC KRAMERS EQUATION 

We consider a system of n interacting particles in the 
3-dimensional space, defined by the Hamiltonian 



n 



P 

2m 



(1) 



where the vectors q = (gi, . . . , (f„) and p = (pi, . . . ,p„) 
indicate the positions and momenta associated to the 
particles, and V{q) is the interaction potential. To sim- 
plify the notation, we assign to each particle the same 
mass m. The dynamics of the system coupled to a heat 
bath at constant temperature T is described by means of 
a Langevin equation 



q = p/m 

p = -\/V + ^2mjTT] - 7p , 



(2) 



where we have fixed the Boltzmann constant fee = 1, the 
friction coefficient is 7, and 77 is a Gaussian white noise: 



M)) = 



(3) 
(4) 



The indices /i and v run over all the configuration space 
degrees of freedom 1, . . . , A^, with N = 3n. 

The phase space probability density VF(q, p, t) evolves 
according to the Kramers equation [fj] 



-i/KT^(q,P,t), 



(5) 



where 
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The Kramers equation can be rewritten as a continuity 
equation for the probability current 



(7) 
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dpij 



dV 



W^(q,p,i). 
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It has been shown [ll|, |13| that a hidden supersymmetry 
is associated with the Kramers equation: By extending 
the space with AN fermion operators 

{af,,al} ^ 6^,1, {bf„bl} = Sf,^ , (8) 

a supersymmetric extension of Eq. ([5]) is obtained 

jjL^U—l ^ fJ,— l 

(9) 

By defining a 2A'^-component vector x such that a;^ ~ 
and XTv+^i = Pfj,, for /i = 1, . . . , the evohition operator 
Eq. ([9]) can be expressed in compact notation 

2N 

HsK ^Hk+J2 ^^jclc, , (10) 

where (ci, . . . , C2jv) = (ai, . . . , ajv, m6i, . . . , m^Af), and 
the matrix A is 

(0 —Snu/m \ 

The solution to the supersymmetric version of Eq. ^ 
can be expressed in the form 

|^W(x,i))= E Vn,...,.Jx,t)c,V..clj-), (12) 

where the function (x, i) has the physical mean- 

ing of probabihty density in the phase space, |— ) is the 
fermion vacuum, and k is the fermion number, that is, 
an eigenvalue of the operator A^f = X^i^Jci. By using 
this notation, the supersymmetric extension of Eq. ([5]) is 
written as 

|^|V('=)(x,i)) =-i/sK|V^'=)(x,t)). (13) 



II. SUPERSYMMETRIC MOLECULAR 
DYNAMICS 

Let us first consider the solution to Eq. in the 
zero-fermion sector, where \^^^'>{x,t)) — W{K,t)\—) . In 
this case we simply recover the Kramers equation ([5l). 
If we start from some initial condition |i/)(x, 0)), we can 
expand the generic state |i/;(x, <)) into right eigenvectors 
of the operator H}^ 

|^(x,t))=5^C„(t)|V'^(x)), (14) 



where //"kIV'q) = ^ali'a)- ^ increases, this sum is 
obviously more and more dominated by the eigenvectors 
with the smallest eigenvalues. For t oo, only the sta- 
tionary state (defined by A = 0) survives. If the system 
is characterized by the presence of two (or more) well 
separated time-scales Tfast ^ TsIow, a corresponding gap 
is also present in the spectrum of Hk- 

It follows that at a time t such that Tfast <^ t <ti Tgiow, 
the evolution of the system is well approximated by a 
linear superposition of the K right eigenvectors below 
the gap: 

K-l 

|V'(x,0)« EC„(0)e-^°*|^«(x)). (16) 

Q=0 

In the framework of the master equation formula- 
tion of non-equilibrium statistical mechanics, it can be 
proved [3, [3, [3] that K suitable linear combinations 
of the right eigenvectors below the gap exist such 
that the associated probabihty densities W{c[, p, t) are 
positive normalized distributions, nonzero only on non- 
overlapping regions of the configuration space, and sta- 
tionary on time-scales much shorter than Tgiow One can 
therefore use these states for a rigorous and general def- 
inition of metastability. It is important to stress that 
this results hold true independently on the origin of the 
time-scale separation. 

The probability distribution W{x, t) can be used to 
define a dynamic free energy: 

r d^^x 

^{t) = J -^[H{x)W{^,t)-^TW{^,t)\nW{K,t)] . 

(17) 

For i — > cx) the probability distribution tends to the 
Boltzmann distribution 

l\mW{ci,p,t) = lexp (^-^^^) , (18) 

where Z is the partition function. It follows that 

lim J^(t) = -rinZ, (19) 

t — >co 

that is the equilibrium definition of the Helmholtz free en- 
ergy in the canonical ensemble. Usual constant tempera- 
ture MD simulations are limited to the study of the zero- 
fermion sector; for a given potential V{q) the Langevin 
equation is numerically integrated for a time t 3> Tgiow 
large enough to reach equilibrium. 

In the one-fermion sector, on the other hand, the wave- 
function Eq. p2)l reads 

27V 

|^W(x,t))=E^,(x,t)ct (20) 
1=1 

and Eq. may be written as 



a 

SO that Eq. (O yields 

|V^(x,t))=5^a(0)e-^°*|^«(x)), (15) 

a 



— V'»(x,i) = -ilK^»(x,i) -E^»,Vj(x,t) (21) 
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where we have made use of Eq. ([10)) . This equation can 
be solved with the ansatz 'ipi{x,t) — ip{'x.,t)wi(t), where 
vir is a vector of dimension 2N that does not depend on 
X = (q, p), and ip(x, t) evolves with the Kramers equation 



d_ 

dt 



(22) 



This leaves for the vector w the evolution equation 

27V 



dt 



(23) 



In order to avoid a divergence of the norm of w, Eq. (123 
can be modified by adding a term: 



= Af{w)wi - ^ Aij wj . 



2N 



(24) 



The norm |w| is now constant provided that we choose 



AA(w) = 



(25) 



The joint distribution T4^(x, w,t) evolves according to 
dW 



dt 



= [-Hk - AA(w)+ 



2N „ / 2N 
O 




- A/'(w)i 



W . 



(26) 



as can be checked by defining 

i;i{x,t) ^ J d^'^^vw^W{x,w,t) (27) 

and integrating by parts. 

The rules of SuSy MD are easily read from the RHS 
of Eq. ([26]) . We are going to explore the free energy 
landscape by means of "walkers" moving around in the 
phase space according to the usual Langevin dynamics 
(first term) and each walker carries a "compass" w which 
evolves with Eq. ([M|) (third term). The second term tells 
us that the number of walkers grows or decreases with 
rate — A/'(w). 

How does the presence of different time-scales reflect 
in the 1-fermion sector of the spectrum of iJsK? In a 
simplified setting where entropy plays no role and the 
separation of time-scales is purely due to the character- 
istics of the energy landscape, the use of a WKB tech- 
nique in the limit T — > shows explicitly [lo| that, while 
the 0-fcrmion states are Gaussians centered on the local 
minima of the energy, the correspondent (i.e. related by 
the supersymmetry) 1-fermion states are the "reduced 
current" densities [11| (obtained by applying the SuSy 
charge operator to the probability currents concen- 
trated on the saddles that separate those minima. In 
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FIG. 1: The free energy profile at the folding temperature, as 
obtained from equilibrium MD simulations. Representative 
structures are shown for the helix-shaped native state and 
the unfolded state. Each contour marks an increase of free 
energy of 1 Kcal/mol. 



other words, the dynamics given by Eqs. ^ [24]) evolves 
in such a way that the walkers quickly (that is, on a 
time-scale larger than rfast but much smaller than Tgiow) 
organize themselves into trails going from one local min- 
imum to another one by overcoming the energy barrier 
along the reaction path Q . 

Since in the zero-fermion sector the right eigenvectors 
below the gap define the metastable states independently 
on the physical source of metastability, it is tempting to 
speculate that the interpretation of 1-fermion low-lying 
states as reaction paths holds also for the general case 
involving entropy, with a single reaction path in the free 
energy landscape standing now for a collection of paths 
in the phase space. As a matter of fact, in the zero- 
temperature limit the dynamic free energy Eq. (|17p re- 
duces to the energy, therefore we can think of the WKB 
argument in Ref. [ll| as a rigorous proof, albeit given in 
a limiting case, of a more general statement. While the 
generalization of the proof to finite temperatures is cur- 
rently in progress, we support here the validity of these 
ideas by showing that indeed SuSy MD can be used to 
efficiently identify reaction paths and saddle points on a 
free energy landscape, in a system where both entropic 
and energetic factors play a role. 



III. THE HELIX-COIL TRANSITION 

The choice of the helix-coil transition as test system 
is a natural one: It is a simple phenomenon, theoreti- 
cally well understood [13], whose free energy is shaped 
by the competition between energy and entropy into a 
landscape with two well defined minima, corresponding 
to the folded and unfolded states (see Fig. [T]). At the 
transition temperature T{ the two minima are equally 



FIG. 2: Relative walkers' density ai t = for two very differ- 
ent initial conditions, that lead to the same final result. The 
two initial distributions of walkers have been generated by: 
(Case A, shown in the left panel) Performing long equilib- 
rium simulations around the transition temperature Tf, and 
(Case B, shown in the right panel) performing very short un- 
folding simulations at a much higher temperature T S> Tf. 



populated. By using a coarse-grained ofF-lattice model 
we keep relatively low the dimensionality of the phase 
space (72 degrees of freedom for our 12-monomer chain, 
see the Appendix for detail), thus reducing the computa- 
tional effort while retaining the relevant physical features 
of a typical two-state folder. 

In order to visualize the results, we need to project 
the 72-dimensional phase space associated to our model 
onto a lower dimensional space spanned by a few reaction 
coordinates ^. Although the definition of appropriate 
reaction coordinates for the characterization of multidi- 
mensional biophysical processes is in general an area of 
active research [18| , a fairly natural set of coordinates is 
associated to the simple helix-coil transition considered 
here: The root mean square deviation (rmsd) [20| 
from the native state x(°) 



rmsd(x,x(")) = min i|(i?x - x(o))| 



i?eso(3) 



(28) 



and the "helicity" 



Nn 



E 



0f))V(A^R-3) 



(29) 



where N^i is the number of residues and , 4>^f'^ are the 
dihedral angles of a generic configuration and of the na- 
tive configuration, respectively. 

As the position x is projected onto the space spanned 
by the reaction coordinates ^, so is the vector w, by 
means of the Jacobian matrix: 



(30) 



where u) is the projected vector. 

The study of the system by means of SuSy MD first 
requires the generation of an initial distribution of a large 
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FIG. 3: The walkers' distribution obtained by SuSy MD (as 
describd in the text) identifies the transition state region and 
reaction paths for the helix-coil transition. The red arrows 
illustrate the orientation of the compasses associated to the 
walkers, and are superimposed to the independently deter- 
mined free energy profile, at the transition temperature Tf. 
Each contour marks an increase of free energy of 1 Kcal/mol. 
A logarithmic scale has been adopted to improve the read- 
ability of the figure: If an arrow in the picture is twice longer, 
the actual norm of the vector is ten times larger. 



number of walkers in the accessible phase space. Each 
walker (x, w) is then evolved independently according to 
Eqs. ^ and Moreover, after each time step 5t, 

there is a probability |A/'(w)(5t| for every walker of being 
eliminated if A/'(w) > or cloned if A/'(w) < 0. 

Figure [2] shows two different distributions of walkers' 
initial configurations that have been used in this study. 
The distributions were generated by running MD simula- 
tions in very different conditions: The initial distribution 
of walkers shown in the left panel (case A, in the follow- 
ing) is obtained by performing equilibrium simulations 
around the folding temperature T ~ Tf, over a very long 
timescale At > TsIow, so that the initial walkers' distri- 
bution mirrors faithfully the free energy landscape. On 
the contrary, the distribution shown in the right panel 
(case B) corresponds to configurations sampled during 
very rapid (At < Tfast) unfolding simulations at a tem- 
perature T ^ Tf. Our experience is that the initial dis- 
tribution of walkers does not affect the result, as long as 
the region of the landscape between the folded and the 
unfolded state is fairly populated. 

In principle, the transition region is simply revealed by 
the alignment of the compasses: They are randomly ori- 
ented within the states, while along the transition path 
they display coherent behavior. In practice, one needs to 
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sift the points according to criteria such as the walkers 
density and the average rate A/'(w) . After thorough test- 
ing, we have selected an analysis protocol consisting of 
the following three steps: 

1. Select a time window; 

2. Select a density threshold; 

3. Select a threshold value for the variance of ■& (de- 
fined below). 

In the following section we detail each step, showing all 
the phases of the process which leads from the raw data 
to the emergence of the reaction path. 

The final result of our analysis is summarized in Fig. [3l 
where the selected walkers are superimposed on the free 
energy profile independently determined by means of ex- 
tensive MD simulations and standard techniques. Re- 
markably, the walker positions and the orientation of 
their compasses clearly highlight the minimum free en- 
ergy path connecting the native and the unfolded states. 
With a more restrictive choice of the various threshold 
values, the transition state region can be pinpointed as 
well. As predicted, the simulation time needed by the 
walkers to find the path is of the order of 10"* time steps, 
significantly shorter than the characteristic time associ- 
ated to activation process, which is around 10^ time steps 
for the helix-coil transition considered here. We envision 
the time separation to be even more pronounced for more 
complex systems. 



IV. DATA ANALYSIS PROCEDURE 

A. Effective migration and time window choice 

The reduced current we want to observe requires a time 
larger than rtast (although much smaller than TsIow) to 
form. On the other hand the current disappears once the 
equilibrium is achieved. A look to the evolution of the 
walker density helps fixing the most profitable time win- 
dow. As an example, we show in Fig. 2] several snapshots 
of the walker distribution obtained starting from the two 
different initial conditions displayed in Fig. [21 While case 
A reflects the Boltzmann distribution at T{, case B is 
quite far from equilibrium. Figure [?] compares the time 
evolution of the walkers' density in the two cases. Fi- 
nally, Fig. [5] shows that when the equilibrium is reached 
any difference due to the different initial conditions is 
lost. 

Based upon the inspection of the walkers' migration, 
we select as time windows the interval [0.5, 6] ps for case 
A, and [0.09, 0.36] ps for case B. The difference in the 
time scales of the walker's migration is due to the fact 
that initial condition B is at higher temperature than 
initial condition A. 




rmad (A) rrnsd (Ji) 



FIG. 4: An effective migration of the walkers is observed on a 
time-scale rfast < t < Taiow, and it is signaled by the changes 
in the relative density. Results shown here correspond to the 
walker density in different time windows, for the initial condi- 
tion A (left figures) and B (right figures), as defined in Figure 
2. 



B. Density threshold and rate distribution 

Once the time window is chosen, in order to reduce the 
unavoidable noise present in the data, we filter out all 
the points of the grid that are not consistently populated 
(i.e. have a density below a given threshold) during the 
migration process. The distribution of the walkers' pop- 
ulation in the space spanned by the reaction coordinates 
shown in Fig. [S]is given by all the walker configurations 
visited during all the independent simulations performed 
within the considered time window. In order to make the 
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FIG. 5: Over long time-scales t ~ TsIow the distribution of 
walkers reach equilibrium and the difference in the initial con- 
ditions is completely lost. The left figure shows the walker 
density at equilibrium for the initial condition A while the 
right figure corresponds to the initial condition B (see Fig.[2|. 
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FIG. 6: Walkers' density averaged over all the time snapshots 
in the chosen time window, as discussed in the text. The 
left figure shows the walker density obtained for the initial 
condition A while the right figure corresponds to the initial 
condition B (see Fig. [2]). 



choice of the density threshold somewhat less arbitrary, 
we adopt the following criterion: the threshold should be 
low enough that we do not disconnect the two metastable 
states, but high enough to have a fair statistics at each 
point of the grid. Within these two boundaries, we veri- 
fied that the actual value of the threshold does not affect 
the final result. After inspection of Fig. [51 we choose the 
value 18 as density threshold for case A and 16 for B. The 
selected configurations cover the native state, the tran- 
sition path and the unfolded state. Now we need some 
quantity to discriminate between the states and the re- 
action path. 

This is a good place to explain the mechanism of 
the walkers' effective migration. If we picture the av- 
erage rate N'{'w) for donation/destruction (defined in 
Eq. ([H])), we notice that the probability of donation is 
larger in the unfolded state region (Fig. [7]). This drives 
the effective migration. One may notice that the rate 
is everywhere negative: In fact, our implementation of 
the supersymmetric Langevin equation is characterized 
by the fact that the number of walkers grows exponen- 
tially. A random decimation of walkers when their num- 
ber exceeds some maximum value is sufficient to solve 
the problem and does not introduce any significant bias 
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FIG. 7: Distribution of the rate A/'(w) for donation or de- 
struction of the walkers, averaged over all the snapshots, on 
the 2-dimensional space spanned by the reaction coordinates. 
Results shown in the left figure are obtained with the initial 
condition A while the right figure corresponds to the initial 
condition B (see Fig. [2]). 



in the final result. 



C. The reaction path revealed 

The walker density alone will not reveal the informa- 
tion we are most interested in: that is, where the re- 
duced current is stronger. This information is stored in 
the "compasses" associated to the walkers: we expect the 
vectors to be strongly coUinear in correspondence of the 
reaction path, and disordered within the states. The av- 
eraged value of the vector is not a reliable quantity to look 
at, because it can be affected by cancellations between 
vectors with same direction but different sign. It is con- 
venient to define the direction angle = arctan(^2/Ci), 
where ^1,^2 are the reaction coordinates we are using: 
Root mean square deviation and helicity, respectively. 
The variance of is a good measure of the coherence be- 
tween the directions of vectors in the same cell of our grid. 
Figure [5] shows that the variance is indeed a good marker 
for the reaction path. By combining the information of 
Figures [SI El [S] we select a set of walkers corresponding 
to densely populated regions, with an associated high 
donation rate, and with a small variance of 

Once a reasonable threshold is chosen for the variance 
of d, the resulting configurations can be compared with 
the free energy profile computed for the same model with 
traditional MD techniques. With a threshold value of 
0.52 for case A and 0.45 for case B, the result is shown in 
Fig. [HI Each vector displayed in the figure at a given po- 
sition (^1,^2) represents an average over a small volume 
d£,id^2 centered in (^i,$2)- AH the figures were obtained 
with the values d^i = 0.04, and d^2 = 0.04. 



V. CONCLUSIONS 

The results presented in this paper show that a su- 
persymmetrically enhanced version of molecular dynam- 
ics can be efficiently used to identify transition states 
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FIG. 8: The average variance of the direction angle ■& — 
arctan(^2/Ci) (where ^1,^2 are the reaction coordinates) indi- 
cates the region where the vectors associated to the walkers 
are more aligned. Results shown in the left figure are obtained 
with the initial condition A while the right figure corresponds 
to the initial condition B (see Fig. [2]). 
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FIG. 9: The compasses of the selected walkers are compared 
with an independently derived free energy profile. The left 
figure shows the results obtained for the initial condition A 
while the right figure corresponds to the initial condition B 
(see Fig. [5)). The red arrows illustrate the orientation of the 
compasses associated to the walkers, and are superimposed 
to the independently determined free energy profile, at the 
transition temperature Tf. Each contour marks an increase 
of free energy of 1 Kcal/mol. Since there are big differences 
in the length of these vectors (as discussed in the text) a 
logarithmic scale has been adopted to improve the readability: 
If an arrow in the picture is twice longer, the actual norm of 
the vector is ten times larger. 

and reaction paths in models of macromolecular sys- 
tems characterized by a clear separation of time-scales 
Tsiow ^ Tfast ■ The great advantage of the method is that 
the simulation does not need to extend over the long 
time-scale Tgiow, since the SuSy Kramers spectrum con- 
tains from the very beginning all the information about 
the topology of the phase space [ll| . The trade-off is that 
instead of a single trajectory, a large number of walkers 
are used to explore the phase space. However, since each 
single walker trajectory is extremely short, SuSy MD is 
easily and efficiently implemented in a parallel computing 



framework. 

Althought the thoretical/mathematical foundation of 
the SuSy MD approach has some similarities with 
the recently prop osed Finite Temperature String (FTS) 
method [22, [2^, there are important differences. A 
comparison of these methods clearly highlights relative 
strengths and weaknesses. While the FTS technique (as 
well as the transition path sampling [24i] that similarly 
relies on evolving a string rather than a point-like object 
in the phase space) requires the definition of an initial 
and a final state, the SuSy walkers are able to find their 
own way without any previous knowledge of the configu- 
rational landscape. In addition, the SuSy approach does 
not require the FTS assumption that the isocommittor 
surface of the reaction could be locally approximated by 
a hyperplane. On the other hand, FTS-based approaches 
bypass the pro blems related to the choice of reaction 
coordinates since they work directly in the high- 
dimensional phase space. Although the work presented 
in this paper was based on the a priori knowledge of a 
good set of reaction coordinates, nothing in the method 
itself require such a step. It should be possible to modify 
the data analysis procedure in such a way that clusters of 
configurations along the reaction path are read directly 
in the phase space. We believe this to be a promising 
direction for further research. 
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APPENDIX: EXPLICIT EXPRESSION OF THE 
POTENTIAL AND VALUES OF THE 
PARAMETERS 

The system selected for our study is a Go-like [2^ 
model that uses a short alpha-helical segment as native 
structure. In particular, we chose the first 12 residues 
of chain A of the Alanine-zipper described in Ref. [2^ 
(PDB code IJCD) as native helical structure. The po- 
tential energy associated to the system is in the form: 
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ATr-I 



Nti-3 



V 



+ei E ' 



5( ^ 



61^ 



1=1 



(A.l) 



r 



where: 



• The residues are numbered from 1 to Nb. and their 
position is represented by the Ca atoms; 

• rij is the distance between residues while = 

• 9i is the angle between the vector from residue i to 
i + I and the vector from i + 1 to i + 2; 

• ipi is the dihedral angle formed by the residues + 
l,i + 2,i + 3; 

• denotes a sum over pairs of residues i,j with 
j - i > 4; 

• C is the native contact map, that is the list of 
residue pairs that are in contact in the native struc- 
ture; 

• The constants Tj-"'' , , ipf^"^ are fixed by the corre- 
sponding values in the native structure; 

• The parameters cTy are set equal to the distance 
between the Ca atoms of residue i and j in the 
native structure. 



The values of the remaining constants in Eq. (jA.ip have 
been chosen as follows: 



k,= 100 kcal/mol/A^ 


kip''— 1 kcal/mol 


£1=5 kcal/mol 


fce=20 kcal/mol/rad^ 


(2) 

k\p —0.5 kcal/mol 


62=1 kcal/mol 



while fjo is equal to 3 A, and the native contact map C 
is 



{(1, 5), (2, 6), (3, 7), (4, 8), (5, 9), (6, 10), (7, 11), (8, 12)} . 



The transition temperature Tf of the system is defined 
as the temperature corresponding to a peak in the heat 
capacity curve. With our choice of the parameters we 
obtain Tf ~ Q.2ei/kB- AH the thermodynamic quanti- 
ties (including the free energy surface reported in Fig. [3]) 
were obtained by combining extensive MD simulations 
at different temperature with the Weighted Histogram 
Analysis Method (WHAM) [13, HI- 

In our implementation, the Langevin equation is solved 
by means of the second-order quasi-symplectic integrator 
described in Ref. [1^ , while the conservation of the norm 
of the vector w is achieved by applying the implicit mid- 
point rule (see, for instance, Ref. 30]). 

The friction coefhcient entering the Langevin equation 
is set to 7 = 2.5 ps~^, and the mass of each particle 
is TO = 100 amu. The time step used in all dynamical 
simulations is 6t ~ 10^'* ps. The SuSy MD simulations 
used 60, 000 independent walkers for each temporal 
snapshot considered. 
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